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Abstract 



Sparse elimination exploits the structure of polynomials by measuring their complexity in terms of 
Newton polytopes instead of total degree. The sparse, or Newton, resultant generalizes the classical 
CZ2 ■ homogeneous resultant and its degree is a function of the mixed volumes of the Newton polytopes. 

We sketch the sparse resultant constructions of Canny and Emiris and show how they reduce the 
problem of root-finding to an eigenproblem. A novel method for achieving this reduction is presented 
which does not increase the dimension of the problem. Together with an implementation of the sparse 
resultant construction, it provides a general solver for polynomial systems. We discuss the overall 
implementation and illustrate its use by applying it to concrete problems from vision, robotics and 
00 ■ structural biology. The high efficiency and accuracy of the solutions suggest that sparse elimination 

1 may be the method of choice for systems of moderate size. 



^ ! 1 Introduction 



The problem of computing all common zeros of a system of polynomials is of fundamental importance 
in a wide variety of scientific and engineering applications. This article surveys an efficient method 
based on the sparse resultant for computing all isolated solutions of an arbitrary system of n polynomials 
in n unknowns. In particular, we exploit the algorithms of Canny and Emiris [5, 11] for constructing 
sparse resultant formulae which yield nontrivial multiples of the resultant. We show that the matrices 
obtained allow the reduction of the root-finding problem to the eigendecomposition of a square matrix. 
The emphasis here is placed on practical issues and the application of our implementation to concrete 
problems. 

We describe very briefly the main steps in sparse elimination and the construction of sparse resultant 
matrices. Most proofs are omitted but can be found in [5, 8, 11, 9, 10]. The study of coordinate rings of 
varieties in K n , where K is a field, has been shown to be particularly useful in studying systems of poly- 
nomial equations. We concentrate on zero-dimensional varieties for which it is known that the coordinate 
ring forms a finite-dimensional vector space and, actually, an algebra over K. An important algorithmic 
question is the construction of an explicit monomial K-hasis for such a space. Based on monomial bases, 
we may generate generic endomorphisms or multiplication maps for any given polynomial, as outlined in 
section 2. 



*Most of this work was conducted as part of the author's Ph.D. thesis in the Computer Science Division of U.C. Berkeley 
(completed in 1994). 



Root finding is reduced to an eigenproblem and then existing techniques are employed from numerical 
linear algebra. An important feature of our method is precisely that it reduces to matrix operations for 
which relatively powerful and accurate implementations already exist. Section 3 discusses this method 
in connection to both resultant algorithms for the case of adding an extra u-polynomial to obtain an 
overconstrained system. This is the classical method, used in defining the u-resultant; it possesses the 
advantage that the matrices have a lot of known structure. 

A relatively novel approach that keeps the number of polynomials fixed is proposed in section 4 where 
one of the variables is hidden in the coefficient field, thus producing an overconstrained system. We show 
how the calculation of all isolated roots again reduces to an eigenproblem. This technique keeps the 
number of polynomials fixed, which has been observed to be important in practice. On the other hand, 
it leads to arbitrary matrix polynomials for which we have to calculate all eigenvalues and eigenvectors. 

This approach has been implemented in C by the author and provides, together with an eigenvalue 
solver, a self-contained and fast polynomial solver. Section 5 describes our implementation and dis- 
cusses the practical issues that arise thereof. In particular, we consider issues of numerical stability and 
conditioning of the matrices. 

Our techniques find their natural application in problems arising in a variety of fields, including prob- 
lems expressed in terms of geometric and kinematic constraints. As an empirical observation, polynomial 
systems encountered in robot and molecular kinematics, motion and aspect ratio calculation in vision 
and geometric modeling are characterized by small mixed volume compared to their Bezout bound. The 
complexity of our methods depends directly on this sparse structure, in contrast to Grobner bases. Sparse 
homotopies were proposed in order to exploit the same structure [18, 28], yet they still suffer from accu- 
racy problems and the possibility that some solutions may not be found. Lastly, resultant-based methods 
include a large fraction of offline computation: the first phase, including the construction of the matrix, 
has to be executed only once for a given set of supports. For every specific instance, the coefficients are 
specialized and then the online phase has to be executed. 

We describe in detail two problems from vision, robot kinematics and structural biology. The first 
problem, analyzed and solved in section 6, is a standard problem from photogrammetry. Given infor- 
mation about a static scene seen by two positions of a camera, the camera motion must be computed. 
When the minimum amount of information is available so that the problem is solvable, an optimal sparse 
resultant matrix is constructed and the numerical answers are computed efficiently and accurately by 
our implementation. Our method exhibits competitive speed, as compared to previous approaches, and 
better accuracy. 

The second application, in section 7, comes from computational biology and reduces to an inverse 
kinematics problem. The symmetry of the molecule at hand explains the high multiplicity of the common 
zeros, which leads us to compare the two approaches of defining an overconstrained system, either by 
adding a u-polynomial or by hiding one of the input variables. Both methods are used for different 
instances in order for all roots to be calculated accurately. 

We conclude with some open questions in section 8. 

2 Sparse Elimination 

Sparse elimination generalizes several results of classical elimination theory on multivariate polynomial 
systems of arbitrary degree by considering the structure of the given polynomials, namely their Newton 
polytopes. This leads to stronger algebraic and combinatorial results in general. Assume that the number 
of variables is n; roots in (C*) ra are called toric. By concentrating on C* we may, consequently, extend our 
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scope to Laurent polynomials. We use x e to denote the monomial x^ 1 ■ ■ ■ x^ n , where e = (ei, . . . , e n ) € Z n 
is an exponent vector or, equivalently, an integer lattice point, and n G Z>i. Let the input Laurent 
polynomials be 

/i,...,/ n € ^[aJi^f 1 ,...,^,^" 1 ] = K^aT 1 ] (1) 

where -fT is a field. 

Let .Aj = supp(/j) = {oji , . . . , aim } c ^™ denote the set, with cardinality ^j, of exponent vectors 
corresponding to monomials in fi with nonzero coefficients. This set is the support of ff. 

fi = ^ ' C{jX lJ , Cjj 7^ 0. 

Definition 2.1 T/ie Newton polytope of fi is the convex hull of support A%, denoted Qi = Couv(Ai) C 
R n . 

For arbitrary sets in R n there is a natural associative and commutative addition operation called 
Minkowski addition. 

Definition 2.2 The Minkowski sum A + B of sets A and B in R n is 

A + B = {a + b\ a G A,b e B} cf. 

If A and B are convex polytopes then A + B is a convex polytope. 

Let Vol(A) denote the Lebesgue measure of A in n-dimensional euclidean space, for polytope A C K n . 

Definition 2.3 Given convex polytopes A\, . . . , A n C R ra , there is a unique real-valued function MV(Ai, . . . , A n ), 
called the mixed volume of A±, . . . , A n which has the following two properties. First, it is multilinear with 
respect to Minkowski addition and scalar multiplication i.e. for (i,p £ R>o and convex polytope A' k C W l 

MV(A 1 ,...,iiA k + P A' k ,...,A n ) = l iMV'(A 1 ,...,A k ,...,A n ) + pMV(A 1 ,...,A! k ,...,A n ). 

Second, 

MV(A 1 , . . . , An) = n\ Vol(Ai), when A 1 = ■ ■ ■ = A n . 

Notationally, we use 

MV(Qi, ...,Q n ) = MV(Ai, ■■■,A n ) = MV(h, ...,/„). 

We are now ready to state Bernstein's theorem [2], the cornerstone of sparse elimination, generalized to 
arbitrary varieties. 

Theorem 2.4 [13, sect. 5.5] Given are polynomials f±,...,f n £ K\ x, x ] with Newton polytopes Qi, • • • , Q n . 
For any isolated common zero a G (C*) n , let i(a) denote the intersection multiplicity at this point. Then 
J2ai( a ) ^ ^W(Qi) • • • iQn), where the sum ranges over all isolated roots. Equality holds when all coeffi- 
cients are generic. 

Canny and Rojas have substantially weakened the requirements for equality [6]. A recent result 
extends the bound on non-toric roots. 
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Theorem 2.5 [19] For polynomials fi, ■ ■ ■ , f n £ x x ] with supports A\, ■ ■ ■ , A n the number of com- 
mon isolated zeros in C n , counting multiplicities, is upwards bounded by MV(Ai U {0}, . . . ,A n U {0}). 

Bernstein's bound is at most as high as Bezout's bound, which is simply the product of the total 
degrees, and is usually significantly smaller for systems encountered in real-world applications. 

The sparse or Newton resultant provides a necessary and generically sufficient condition for the exis- 
tence of toric roots for a system of n + 1 polynomials in n variables: 

/i,...,/„+i SKfax- 1 ]. (2) 

To define the sparse resultant we regard a polynomial fi as a generic point Cj = (cji, . . . , Cj mi ) in the space 
of all possible polynomials with the given support Ai = supp(/j), where rrii is the number of nonzero 
terms. It is natural to identify scalar multiples, so the space of all such polynomials contracts to the 
projective space P^ i_1 or, simply, P™ 1-1 . Then the input system (2) can be thought of as a point 

c = (ci, . . . , c n+ i) G P™ 1 " 1 x • • • x p^+i- 1 . 

Let Zq = Zq(A\, . . . ,A n +i) be the set of all points c such that the system has a solution in (C*) n and let 
Z = Z(Ai, . ■ . ,An+i) denote the Zariski closure of Zq in the product of projective spaces. It is proven 
in [25] that Z is an irreducible variety. 

Definition 2.6 The sparse resultant R = R(Ai, ■ ■ ■ , A n +i) of system (2) is a polynomial in Z[c]. If 
codim(Z) = 1 then R(A\, . . . , An+i) is the defining irreducible polynomial of the hypersurface Z. If 
codim(Z) > 1 then R(Ai, . . . ,A n +i) = 1. 

Let degy. R denote the degree of the resultant R in the coefficients of polynomial fi and let 

MV-i = MV(Q U Qi-!,Q i+1 , . . . , Q n+1 ) for i = 1, . . . , n + 1. 

A consequence of Bernstein's theorem is 

Theorem 2.7 [25] The sparse resultant is separately homogeneous in the coefficients Ci of each fi and its 
degree in these coefficients equals the mixed volume of the other n Nwton polytopes i.e. degj. R = MV-i. 

Canny and Emiris [5, 8, 11] have proposed the first two efficient algorithms for constructing resultant 
matrices i.e., matrices in the coefficients whose determinant is a nontrivial multiple of the sparse resultant. 
The first algorithm relies on a mixed subdivision of the Minkowski Sum, while the second constructs the 
matrix in an incremental fashion. For those cases where it is provably possible, the incremental algorithm 
yields optimal matrices, so that the determinant equals the resultant. For general systems, it typically 
produces matrices that are at most 3 times larger than optimal. Let 

Q = Qi + — h Qn+i c R n 
be the Minkowski sum of all Newton polytopes. Let the resultant matrix be M. Now let 

£ = (Q + ) n z n 

be the set that indexes the rows and columns of M in a bijective way, where 8 6 Q n is an arbitrarily small 
and sufficiently generic vector. Clearly, M has dimension \£\. The incremental algorithm also indexes 
the rows and columns with monomials, or equivalently, lattice points in £. This algorithm, though, is 
different and may select some points more than once. In practice, this algorithm produces significantly 
smaller matrices but we have no formal result on their dimension. Irrespective of the algorithm applied, 
the resultant matrix M has the following properties. 
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Theorem 2.8 [10] Matrix M obtained by either algorithm is well-defined, square, generically nonsingular 
and its determinant is divisible by the sparse resultant R. 

To be more precise, the rows of M are indexed by a pair composed of a monomial and an input 
polynomial. The entries of the respective row are coefficients of this polynomial. The degree of det M in 
the coefficients of /j for i = l,...,n + lis greater or equal to MV^i. 

To solve system (1) we define an overconstrained system by one of the two ways below. We apply the 
resultant matrix construction on the new system and use the following properties. Let X = I(fi, ... , f n ) 
be the ideal generated by polynomials (1) and V = V(f±, . . . , f n ) € (K ) n their variety, where K is the 
algebraic closure of field K. Generically, V has dimension zero. Then, its coordinate ring if[x,x _1 ]/X is 
an m-dimensional vector space over K by theorem 2.4, where 

m = MV(f 1 ,...,f n ) = MV(Q 1 ,...,Q n ). 

Using the subdivision-based construction it is easy to show [9] that generically a monomial basis of 
K[x, x _1 ]/X can be found among the monomials of 

Qi + ■ • • + Qn C R n . 

Moreover, resultant matrix M produces the multiplication map for any given fo. This is a matrix, or an 
endomorphism, that serves in computing in the coordinate ring and essentially allows computation of the 
common roots of /i = ••• = /„ = 0. 

3 Adding a Polynomial 

The problem addressed here is to find all isolated roots a G V where V C (K*) n is the zero-dimensional 
variety of (1), with cardinality bounded by m = MV(Q\, . . . , Q n ). In addition to zero-dimensional, the 
ideal X = X(/i, . . . , /„) is assumed to be radical, or self-radical i.e., X = y/T, which is equivalent to saying 
that all roots in V are distinct. This requirement is weakened later. 

An overconstrained system is obtained by adding extra polynomial /o to the given system. We choose 
fo to be linear with coefficients coj and constant term equal to indeterminate u. 

/ = n + c iXiH hc „x n G 

Coefficients coj, j = 1, . . . , n, should define an injective function 

fo : V -> K : f Q (a). 

There are standard deterministic strategies for selecting coj so that they ensure the injective property. 
In practice, coefficients coj may be randomly distributed in some range of integer values of size S > 1, 
and a bad choice for coi , . . . , co n is one that will result in the same value of fo — u at two distinct roots 
a and a'. Assume that a and a' differ in their i-th coordinate for some i > 0, then fix all choices of con- 
fer j ^ i; the probability of a bad choice for co% is 1/5", and since there are (™) pairs of roots, the total 
probability of failure for this scheme is 

Prob[failure] < \\ / S : c 0j G {1, . . . , S}, j = 1, . . . , n. 
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It suffices, therefore, to pick coj from a sufficiently large range in order to make the probability of success 
arbitrarily high. Moreover, it is clear that any choice of coefficients can be tested deterministically at the 
end of the algorithm. 

Either algorithm for the resultant matrix may be used to build matrix M. As before, the vanishing 
of det M is a necessary condition for the overconstrained system to have common roots. For a G V, u is 
constrained to a specific value determined by the coj coefficients. The construction of M is not affected 
by this definition of /o- Let monomial set £ index the columns of M and partition M so that the lower 
right square submatrix M22 depends on u and has size r. Suppose for now that the upper left square 
submatrix Mq is nonsingular. Clearly r > m and equality holds if M is obtained by the subdivision 
resultant algorithm. By the construction of M and for an arbitrary a G (K ) n 



Mn M12 
M 2 i M 22 {u) 



q,p G G {0, 1, . . . ,n}, 



(3) 



Now M'(u) = M22(u) — M2iM 11 1 Mi2 where its diagonal entries are linear polynomials in u and 

a pi f (a,u) 



M'{u) 



a Pr f (a,u) 



(4) 



where B = {b\, . . . ,b r } and {p±, . . . ,p r } index the columns and rows, respectively, of M22 and thus M'. 
For a root a G V and for 



3=1 

the right hand side vector in (4) is null. Let v' a = [a bl , . . . ,a br ] and write M'(u) = M' + ul, where M' 
now is numeric and / is the r x r identity matrix. Then 



(M' + ul)v' a = 



M'- [^cojOiAl 



0. 



This essentially reduces root-finding to an eigenproblem since, for every solution of the original system, 
there is an eigenvalue and eigenvector of M and hence of M'. Below we study how to compute a candidate 
solution from every eigenvalue-eigenvector pair. 

If the generated ideal I is radical then every eigenvalue has algebraic multiplicity one with probability 
greater or equal to 1 — (™)/S. We can weaken the condition that I be radical by requiring only that each 
eigenvalue has geometric multiplicity one. This equals the dimension of the eigenspace associated with an 
eigenvalue. If there exist eigenvalues of higher geometric multiplicity this technique fails: then we may 
use the fact that specializations of the n-resultant yield the root coordinates. Alternatively we can define 
an overconstrained system by hiding a variable as in the next section and derive the root coordinates one 
by one. 

In what follows we assume that all eigenvalues have unit geometric multiplicity. Hence it is guaranteed 
that among the eigenvectors of M' we shall find the vectors v' a for a G V. By construction of M [9] each 
eigenvector v' a of M' contains the values of monomials B at some common root a G (K ) n . By (3) we 
can define vector v a as follows: 



M n v a + M 12 v' a = 0^v a = -M^M 12 v' a . 



(5) 
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The size of v a is \£\ — r, indexed by £ \ B. It follows that vectors v a and v' a together contain the values 
of every monomial in £ at some root a. 

Theorem 3.1 [10] Assume £ spans Z™. Then there exists a polynomial-time algorithm that finds a subset 
of n + l affinely independent points in £. Given v a , v' a and these points, we can compute the coordinates 
of root a G V(I). If all n + l independent points are in B then v' a suffices. 

In practice, most of the operations described here are not implemented with (exact) rational arith- 
metic but are instead carried out over floating point numbers of fixed size. An important aspect of this 
computation is numerical error, which we discuss below in a separate section and in the particular context 
of specific applications later. 

Theorem 3.2 [10] Suppose that system (1) generates a zero- dimensional radical ideal, matrix M has been 
computed such that M\\ is nonsingular and n + l affinely independent points in £ have been computed. Let 
r be the size of M' , fi the maximum number of monomials in any support and d the maximum polynomial 
degree in a single variable. Then all common zeros of the polynomial system are approximated in time 

0(|£| 3 + rnVogd). 

It is clear that for most systems the arithmetic complexity is dominated by the first term, namely the 
complexity of matrix operations and in particular the eigendecomposition. Under reasonable assumptions 
the complexity becomes [10] 

2°^m 3 , where m = MV(fi, . . . , f n ). 

As expected, the complexity is single exponential in n and polynomial in the number of roots. 

One hypothesis concerned the nonsingularity of M\\. When it is singular the resultant matrix is 
regarded as a linear matrix polynomial in u and the values and vectors of interest are its singular values 
and right kernel vectors. Finding these for an arbitrary matrix polynomial is discussed in the next section. 



4 Hiding a Variable 

An alternative way to obtain an overconstrained system from a well-constrained one is by hiding one of 
the original variables in the coefficient field. Hiding a variable instead of adding an extra polynomial 
possesses the advantage of keeping the number of polynomials constant. Our experience in solving 
polynomial systems in robotics and vision suggests that this usually leads to smaller eigenproblems. The 
resultant matrix is regarded as a matrix polynomial in the hidden variable and finding its singular values 
and kernel vectors generalizes the u-polynomial construction of the previous section. 
Again we suppose the ideal is zero-dimensional and radical. Formally, given 

f ,...,f n e K[x 1 ,Xi 1 ,...,x n+1 ,x~l 1 ] (6) 

we can view this system as 

/o, •••,/« € K[x„ + i][xi,x ] ; 1 ,...,x n ,x~ 1 ], (7) 

which is a system of n + 1 Laurent polynomials in variables x\, . . . ,x n . Notice that we have multiplied 
all polynomials by sufficiently high powers of x ra +i in order to avoid dealing with denominators in the 
hidden variable x n+ \. This does not affect the system's roots in (K ) n+1 . 

The sparse resultant of this system is a univariate polynomial in x n +i- We show below how this 
formulation reduces the solution of the original well-constrained system to an eigenproblem. 
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Theorem 4.1 Assume that M is a sparse resultant matrix for (7), with the polynomial coefficients 
specialized. Let a = (a±, . . . ,a n ) G (K ) n such that (a, a n +i) G (K ) n+1 is a solution of fi = •■■ = 
fn+i = 0. Then M(a n+ i) is singular and column vector w = [a qi , . . . , a Qc ] lies in the right kernel of 
M{a n+ i), where £ = {q±, . . . ,q c } C Z n are the exponent vectors indexing the columns of M. 

Proof For specialized polynomial coefficients, M(a n +i) is singular by definition. By construction, right 
multiplication by a vector of the column monomials specialized at a point produces a vector of the values 
of the row polynomials at this point. Let the i-th row of M contain the coefficients of x Pj fj, then 



M(a n+ i)w = M(a n+ i) 



' a qi ' 




' a Pl f il (a,a n+1 ) ' 




' " 




. a qc _ 




. a Pc f ic (a,a n+1 ) _ 




. . 


: h 



: h, . . . ,i c £ {0, 1, . . . ,n}. 



□ 



Computationally it is preferable to have to deal with as small a matrix as possible. To this end we 
partition M into four blocks so that the upper left submatrix Mn is square, nonsingular and independent 
of x n+ i. Row and column permutations do not affect the matrix properties so we apply them to obtain 
a maximal Mn. 

Gaussian elimination of the leftmost set of columns is now possible and expressed as matrix multipli- 
cation, where / is the identity matrix of appropriate size: 







-l 



-M 2 i(x„+i)M n 



M u M 12 (x n+1 ) 
M 2 i(x n+1 ) M 22 (x n+ i) 



Mn M 12 (x n+1 ) 
M'(x n+1 ) 



(8) 



where 



M'{x n+l ) = M 22 (x n+1 ) - M 2 i(x n+1 )M 11 1 M 12 (x n+1 ). 



Let B C £ index M' . We do not have a priori knowledge of the sizes of Mn and M' whether the 
subdivision or the incremental algorithm has constructed M. 

Corollary 4.2 Let a = (a±, . . . , a n ) G (K*) n such that (a,a n +i) G (K*) n+l is a common zero of 
fo = ■ ■ ■ = f n = 0. Then det M'(a n+ i) = and, for any vector ?/ = [■■■ ofl ■ ■ ■], where q ranges over B, 

m'k+iV = o. 

To recover the root coordinates, B must affinely span Z n , otherwise we have to compute the kernel 
vector of matrix M which equals the concatenation of vectors v and v', where v' is the kernel vector of 
M' and v is specified from (8): 



Mn Mi 2 (a„+i) 
M'{a n+1 ) 





V 




' " 




v' 








M u v + M 12 (a n+1 )v' = 
^ v = -Mf 1 1 Mi 2 (a n+ i)t; / , 



since Mn is defined to be the maximal nonsingular submatrix. The concatenation [v, v'\ is indexed by 
£ which always includes an affinely independent subset unless all ra-fold mixed volumes are zero and no 
roots exist. Then, we recover all root coordinates by taking ratios of the vector entries. 
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We now concentrate on matrix polynomials and their companion matrices; for definitions and basic 
results consult [15]. Denote the hidden variable by x, then 



fi £ K[x][x 1 ,x 1 ,...,x n ,x n ], i = 0,...,n, 



(9) 



and we denote the univariate matrix M'(x n+ i) by A{x). Let r be the size of A, and d > 1 the highest 
degree of x in any entry. We wish to find all values for x at which matrix 



A(x) 



x d A d + Ai-i + • • • + xAi + A 



becomes singular, where matrices Ad, . . . ,Aq are all square, of order r and have numerical entries. We 
refer to A{x) as a matrix polynomial with degree d and matrix coefficients A{. The values of x that make 
A(x) singular are its eigenvalues. For every eigenvalue I, there is a basis of the kernel of A(l) defined 
by the right eigenvectors of the matrix polynomial associated to I. This is the eigenproblem for matrix 
polynomials, a classic problem in linear algebra. 

If Ad is nonsingular then the eigenvalues and right eigenvectors of A(x) are the eigenvalues and right 
eigenvectors of a monic matrix polynomial. Notice that this is always the case with the it-resultant 
formulation in the previous section. 

A^A{x) = x d I + x^AfAd-! + ■■■ + xA~ d 1 A l + AfA , 

where / is the mx m identity matrix. The companion matrix of this monic matrix polynomial is defined 
to be a square matrix C of order rd. 



C = 






-Aj 1 A, 



o 



-Al l A, 







-Al'Ad-! 



It is known that the eigenvalues of C are precisely the eigenvalues of the monic polynomial, whereas 
its right eigenvectors contain as subvectors the right eigenvectors of A'7 1 A(x). Formally, assume w = 

[v± , . . . , Vd] £ K md is a (nonzero) right eigenvector of C, where each Vi G K m , i = 1, . . . , d. Then v\ is a 
(nonzero) right eigenvector of A^ 1 A(x) and vi = l l ~ 1 v\, for i = 2, . . . , d, where I is the eigenvalue of C 
corresponding to w. 

We now address the question of a singular leading matrix in a non-monic polynomial. The following 
transformation is also used in the implementation in order to improve the conditioning of the leading 
matrix. 



Lemma 4.3 Assume matrix polynomial A(x) is not identically singular for all x and let d be the high- 
est degree in x of any entry. Then there exists a transformation x i-> (t\y + ^2) / (^32/ + £4) for some 
ti,t2,t3,t4 £ Z, that produces a new matrix polynomial B(y) of the same degree and with matrices of the 
same dimension, such that B(y) has a nonsingular leading coefficient matrix. 

It is easy to see that the resulting polynomial B(y) has matrix coefficients of the same rank, for 
sufficiently generic scalars t\,t2,ts,t4, since every matrix is the sum of d+ 1 scalar products of Aj. Thus 
this transformation is often referred to as rank balancing. 
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Theorem 4.4 If the values of the hidden variable in the solutions of fi = ■ ■ ■ = f n +i = 0, which 
correspond to eigenvalues of the matrix polynomial, are associated to eigenspaces of unit dimension and 
A(x) is not identically singular then we can reduce root-finding to an eigenproblem and some evaluations 
of the input polynomials at candidate roots. 

Proof We have seen that the new matrix polynomial B(y) has a nonsingular leading coefficient. More- 
over, finding its eigenvalues and eigenvectors is reduced to an eigenproblem of the companion matrix C. 
By hypothesis, the eigenspaces of C are one-dimensional therefore for every root there is an eigenvector 
that yields n coordinates of the root by theorem 3.1. The associated eigenvalue may be either simple or 
multiple and yields the value of the hidden variable at the same root. Notice that extraneous eigenvectors 
and eigenvalues may have to be rejected by direct evaluation of the input polynomials at the candidate 
roots and a zero test. The right eigenvectors of B(y) are identical to those of A(x) but any eigenvalue I 
of the former yields (t\l + t2)/{t^l + t±) as an eigenvalue of A(x). □ 

The condition that all eigenspaces are unit-dimensional is equivalent to the solution coordinate at 
the hidden variable having unit geometric multiplicity. For this it suffices that the algebraic multiplicity 
of these solutions be one i.e., all hidden coordinates must be distinct. Since there is no restriction in 
picking which variable to hide, it is enough that one out of the original n + 1 variables have unit geometric 
multiplicity. If none can be found, we can specialize the hidden variable to each of the eigenvalues and 
solve every one of the resulting subsystems. 

Our complexity bounds shall occasionally ignore the logarithmic terms; this is expressed by the use 
of 0*(-). Let n the maximum number of monomials in any (n-variate) polynomial of (7) and d the 
maximum polynomial degree in a single variable; this is typically larger than the highest degree of the 
hidden variable but in the worst case they are equal. 

Theorem 4.5 [10] Suppose that the ideal of (7) is zero- dimensional, the coordinates of the hidden variable 
are distinct, matrix M has been computed such that Mu is nonsingular and the resulting matrix polynomial 
A(x) is regular. In addition, a set of affinely independent points in \£\ has been computed. Then all 
common isolated zeros are computed with worst-case complexity 

0*(\£\ 3 d + MM(rd) + rdn 2 fi) = 0*(\£\ 3 d 3 + \£\dn 2 fi). 

Under some reasonable assumptions about the input the complexity becomes [10] 

2°( n '0(m 6 ), where m = MV(f\, f n ). 

It is clear by the first bound on arithmetic complexity that the most expensive part is the eigendecompo- 
sition of the companion matrix. Moreover, the problem of root-finding is exponential in n as expected. 

5 Implementation and Numerical Accuracy 

Our implementation is entirely in Ansi C. The overall method has two stages, one online and one offline. 
A program of independent interest, which makes part of this package, is already available. It is the 
implementation of a fast algorithm for computing mixed volumes [11] and can be obtained by anonymous 
ftp on robotics.eecs.Berkeley.edu, from directory MixedVolume. 

One advantage of the resultant method over previous algebraic as well as numerical methods is that 
the resultant matrix need only be computed once for all systems with the same set of exponents. So this 
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step can often be done offline, while the eigenvalue calculations to solve the system for each coefficient 
specialization are online. 

Another offline operation is to permute rows and columns of the given resultant matrix in order to 
create a maximal square submatrix at the upper left corner which is independent of the hidden variable. 
In order to minimize the computation involving polynomial entries and to reduce the size of the upper left 
submatrix that must be inverted, the program concentrates all constant columns to the left and within 
these columns permutes all zero rows to the bottom. 

To find the eigenvector entries that will allow us to recover the root coordinates it is typically sufficient 
to examine B indexing M22 and search for pairs of entries corresponding to exponent vectors v\ , v 2 such 
that vi — V2 = (0. . . . , 0, 1, 0, . . . , 0). This will let us compute the i-th coordinate if the unit appears at the 
i-th position. For inputs where we need to use points in £ \ B, only the entries indexed by these points 
must be computed in vector v a , in the notation of theorem 3.1. In any case the complexity of this step 
is dominated. 

The input to the online solver is the set of all supports, the associated coefficients and matrix M 
whose rows and columns are indexed by monomial sets. The output is a list of candidate roots and the 
values of each at the given polynomials. The test for singularity of the matrix polynomial is implemented 
by substituting a few random values in the hidden variable and testing whether the resulting numeric 
matrix has full rank. For rational coefficients this is done in modular arithmetic so the answer is exact 
within a finite field, therefore nonsingularity is decided with very high probability. 

Most of the online computation is numeric for reasons of speed; in particular we use double precision 
floating point numbers. We use the LAPACK package [1] because it implements state-of-the-art algo- 
rithms, including block algorithms on appropriate architectures, and provides efficient ways for computing 
condition numbers and error bounds. Moreover it is publicly available, portable, includes a C version 
and is soon to include a parallel version for shared-memory machines. Of course it is always possible to 
use other packages such as EISPACK [26] or LINPACK [3]. 

A crucial question in numerical computation is to predict the extent of roundoff error; see for in- 
stance [16]. To measure this we make use of the standard definition of matrix norm \\A\\ p for matrix A. 
We define four condition numbers for A: 

K P (A) = ll^llpHA -1 !^, p = 1,2, 00, F, A square, and if A singular k p (A) = 00, 

where F denotes the Frobenius norm. K2{A) equals the ratio of the maximum over the minimum singular 
value of A. For a given matrix, the ratio of any two condition numbers with p = 1,2, 00, F is bounded 
above and below by the square of the matrix dimension or its inverse. We also use \\v\\ p for the p-norm 
of vector v. 

As precise condition numbers are sometimes expensive to compute, LAPACK provides approximations 
of them and the associated error bounds. Computing these estimates is very efficient compared to the 
principal computation. These approximations run the risk of underestimating the correct values, though 
this happens extremely seldom in practice. Error bounds are typically a function of the condition number 
and the matrix dimension; a reasonable assumption for the dependence on the latter is to use a linear 
function. Small condition numbers indicate well-behaved or well- conditioned matrices, e.g. orthogonal 
matrices are perfectly conditioned with k = 1. Large condition numbers characterize ill-conditioned 
matrices. 

After permuting rows and columns so that the maximal upper left submatrix is independent of u or 
the hidden variable x, we apply an LU decomposition with column pivoting to the upper left submatrix. 
We wish to decompose the maximal possible submatrix so that it has a reasonable condition number. 
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The maximum magnitude of an acceptable condition number can be controlled by the user; dynamically, 
the decomposition stops when the pivot takes a value smaller than some threshold. 

To compute M' we do not explicitly form M-Q 1 but use its decomposition to solve linear problem 
M\\X = M\2 and then compute M' = M22 — M21X. Different routines are used, depending on re(Mn), 
to solve the linear problem, namely the slower but more accurate dgesvx function is called when this 
condition number is beyond some threshold. Let Xj denote some column of X and let Xj be the respective 
column if no roundoff error were present. Then the error is bounded [16] by 

l|a ~f l|o ° ^ 4 • 10~ 15 (|£:| - r) Koo (M n ), 

Halloo 

where \£\ — r is the size of M\\ and we have used 2 • 1CT 16 as the machine precision under the IEEE 
double precision floating point arithmetic. 

For nonsingular matrix polynomials A(x) we try a few random integer quadruples (^1,^2,^3,^4)- We 
then redefine the matrix polynomial to be the one with lowest n(Ad). This operation of rank balancing 
is indispensable when the leading coefficient is nonsingular as well as ill-conditioned. Empirically we 
have observed that for matrices of dimension larger than 200, at least two or three quadruples should be 
tried since a lower condition number by two or three orders of magnitude is sometimes achieved. The 
asymptotic as well as practical complexity of this stage is dominated by the other stages. 

If we manage to find a matrix polynomial with well-conditioned Ad we compute the equivalent monic 
polynomial and call the standard eigendecomposition routine. There are again two choices in LAPACK 
for solving this problem with an iterative or a direct algorithm, respectively implemented in routines 
hsein and trevc. Experimental evidence points to the former as being faster on problems where the 
matrix size is at least 10 times larger than the mixed volume, since an iterative solver can better exploit 
the fact that we are only interested in real eigenvalues and eigenvectors. 

If Ad is ill-conditioned for all linear tj transformations we build the matrix pencil and call the general- 
ized eigendecomposition routine dgegv to solve Cix + Cq. The latter returns pairs (a, ft) such that matrix 
C\ol + Co/3 is singular. For every (a, P) pair there is a nonzero right generalized eigenvector. For nonzero 
P we obtain the eigenvalues as a/P, while for zero P and nonzero a the eigenvalue tends to infinity and 
depending on the problem at hand we may or may not wish to discard it. The case a = P = occurs if 
and only if the pencil is identically zero within machine precision. 

In recovering the eigenvector of matrix polynomial A(x) from the eigenvector of its companion matrix, 
we can use any subvector of the latter. We choose the topmost subvector when the eigenvalue is smaller 
than the unit, otherwise we use one of the lower subvectors. Nothing changes in the algorithm, since 
ratios of the entries will still yield the root, yet the stability of these ratios is improved. 

There are certain properties of the problem that have not been exploited yet. Typically, we are 
interested only in real solutions. We could concentrate, therefore, on the real eigenvalues and eigenvectors 
and choose the algorithms that can distinguish between them and complex solutions at the earliest stage. 
Moreover, bounds on the root magnitude may lead to substantial savings, as exemplified in [22]. In the 
rest of this article we examine these issues in the light of concrete applications of our program. 

6 Camera Motion from Point Matches 

Camera motion reconstruction, or relative orientation, in its various forms is a basic problem in pho- 
togrammetry, including the computation of the shape of an object from its motion. Formally, we are 
interested in the problem of computing the displacement of a camera between two positions in a static 
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system 


operation 


CPU time 


6x6 


mixed volume 


lm 16s 


6x5 


sparse resultant (offline) 


12s 


6x6 (first) 


root finding (online) 


0.2s 


6x6 (second) 


root finding (online) 


Is (Sun Sparc 20) 



Table 1: Camera motion from point matches: running times are measured on a DEC Alpha 3000 except 
for the second system which is solved on a Sun Sparc 20. 



environment. Given are the coordinates of certain points in the two views under perspective projection 
on calibrated cameras. Equivalently, this problem consists in computing the displacement of a rigid body, 
whose identifiable features include only points, between two snapshots taken by a stationary camera. We 
consider, in particular, the case where the minimum number of 5 point matches is available. In this case 
the algebraic problem reduces to a well-constrained system of polynomial equations and we are able to 
give a closed-form solution. 

Typically, computer vision applications use at least 8 points in order to reduce the number of possible 
solutions to 3 and, for generic configurations, to one. In addition, computing the displacement reduces 
to a linear problem and the effects of noise in the input can be diminished [20]. Our approach shows 
performance comparable to these methods and, as it requires the minimum number of data points, it is 
well-suited for detecting and eliminating outliers in the presence of noise. These are data points which 
are so much affected by noise that they should not be taken into account. 

6.1 Algebraic Formulation 

To formalize, let orthonormal 3x3 matrix R € SO(3,R) denote the rotation. Let (column) vector t 6 M 3 
denote the camera translation in the original frame of reference. The 5 points in the two images are 
(column) vectors a*, a' { € P|, for i = 1, . . . , 5 in the first and second frame of reference respectively. It is 
clear that the magnitude of t cannot be recovered, hence there are 5 unknowns, 3 defining the rotation 
and 2 for the translation. 

The following quaternion formulation was independently suggested by J. Canny and [17]. Let x,y,z 
be 3-vectors and x = [xo,x],y = [yo,y],z = [zq,z] 6 M 4 be arbitrary quaternions. Let xy represent a 
quaternion product and z* = [zq,—z] be the conjugate quaternion of i. Quaternions q = [qo,q], t = 
[to,t] £ R 4 represent the rotation and translation respectively. A rotation represented by angle (j> and 
unit 3- vector s is expressed uniquely by quaternion [cos (j>/2, sin (j>/2 s], hence any rotation quaternion has 
unit 2-norm. 

The equations in terms of the quaternions are homogeneous. After dehomogenization we obtain 6 
polynomials in 6 variables organized in two 3-vectors. We denote these vectors by q, d. 

(a] q){d 7 a-) + a- + (a, x qf a- + (a, x qf \d x a-) + a] \d x a-) = 0, i = l,...,5 

\-d T q = (10) 

This system is well-constrained so we can apply Bernstein's bound to approximate the number of 
roots in (C*) 6 . The mixed volume of this system is 20, which is an exact bound [7]. The performance of 
our implementation on mixed volume is shown in table 1 for the DEC Alpha 3000 of table 2. 
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machine 


clock rate [MHz] 


memory [MB] 


Spcclnt92 


SpecFP92 


DEC Alpha 3000/300 


150 


64 


67 


77 


Sun Sparc 20/61 


60 


32 


95 


93 



Table 2: Hardware specifications 



6.2 Applying the Resultant Solver 

Resultant matrix M is of dimension 60, while only 20 columns contain hidden variable q\. 40 x 40 
submatrix M\\ is inverted and the resulting 20 x 20 pencil is passed to an eigendecomposition routine. 
The running times for the two main phases are reported in table 1. The rest of this section examines the 
accuracy of our procedure on two specific instances from [12]. We use error criterion 

5 1 

E ium .nM KtxRa^, (11) 

i=l H^ll ' H a iH 

where | • | denotes absolute value. || • || is the vector 2-norm and t, R are the calculated translation and 
rotation. This expression vanishes when the solutions to R, t are exact, otherwise it returns the absolute 
value of error normalized by the input norm. 

The first example admits the maximum of 10 pairs of real solutions. M\\ is well-conditioned with 
k < 10 3 , so it is factorized to yield an optimal pencil of dimension 20. The latter has a well-conditioned 
leading coefficient with k < 10 4 , so it yields a monic polynomial on which the standard eigendecomposition 
is applied. The maximum value of error criterion (11) is less than 1.3 • 10~ 6 , which is very satisfactory. 

The input parameters are sufficiently generic in order to lead to the maximum number of real solutions. 
This explains the good conditioning of the matrices. The second example looks at a less generic instance: 
namely a configuration that has only 8 distinct real solutions which are clustered close together. This is 
manifest in the fact that matrix M\\ is ill-conditioned. The input is constructed by applying a rotation 
of 60.00145558° about the z axis in the first frame and a translation of (.01, .01, —1). 

Exactly the same procedure is used as before to produce a 60 x 60 resultant matrix, but now the upper 
leftmost 40 x 40 submatrix M\\ has k > .6- 10 5 so we choose not to invert it. Instead, the 60 x 60 pencil is 
considered: the leading matrix has k > 10 8 in its original form and for any of 4 random transformations, 
hence it is not inverted and the generalized eigendecomposition is applied. 

Applying error criterion (11), the largest absolute value is 7.4 • 10~ 5 , hence all 16 real solutions are 
quite accurate. There are another 4 complex solutions and 40 infinite ones. The total CPU running time 
for the online phase is, on the average, 1 second on the Sun Sparc 20 of table 2. 

It is interesting to observe in connection to roots with zero coordinates, that in the second example 
our solver recovers roots in C n \ (C*) n , namely we find a camera motion whose rotation quaternion q 
has qi = q2 = 0. Such roots can be thought of as limits of roots in (C*) n as the system coefficients 
deform. As long as the variety does not generically reside in C n \ (C*) n , roots with zero coordinates will 
always be recovered. This is typically the case in practical applications. For the particular example, 
there is a stronger reason why all roots are recovered, namely all polynomials include a constant term 
(see theorem 2.5). 

There exist various implementations of linear methods requiring at least 8 points, including Lu- 
ong's [21]. We have been able to experiment with this program which implements the least-squares 
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method of [27], and found it faster on both instances but less accurate on the second one. In particular, 
we chose specific solutions to generate an additional 3 matches for each of the problems above. The 
average CPU time on the two examples is 0.08 seconds on a Sun Sparc 20. On the first example the 
output is accurate to at least 7 digits. On the second example, Luong's implementation returns a rotation 
of 60.0013° about (-10" 5 , 10" 5 , 1) and a unit translation vector of (-10~ 5 , -10~ 5 , -1). But the latter 
differs significantly from the true vector t = (.01, .01, —1). 

7 Conformational Analysis of Cyclic Molecules 

A relatively new branch of computational biology has been emerging as an effort to apply successful 
paradigms and techniques from geometry and robot kinematics to predicting the structure of molecules, 
embedding them in euclidean space and finding the energetically favorable configurations [24, 10]. The 
main premise for this interaction is the observation that various structural requirements on molecules can 
be modeled as geometric or kinematic constraints. 

This section examines the problem of computing all conformations of a cyclic molecule, which reduces 
to an inverse kinematics problem. Conformations specify the 3-dimensional structure of the molecule. 
It has been argued by Go and Scherga [14] that energy minima can be approximated by allowing only 
the dihedral angles to vary, while keeping bond lengths and bond angles fixed. At a first level of ap- 
proximation, therefore, solving for the dihedral angles under the assumption of rigid geometry provides 
information for the energetically favorable configurations. 

We consider molecules of six atoms to illustrate our approach and show that the corresponding 
algebraic formulation conforms to our model of sparseness. Our resultant solver is able to compute all 
solutions accurately even in cases where multiple solutions exist. 

7.1 Algebraic Formulation 

The molecule has a cyclic backbone of 6 atoms, typically of carbon. They determine primary structure, 
the object of our study. Carbon-hydrogen or other bonds outside the backbone are ignored. The bond 
lengths and angles provide the constraints while the six dihedral angles are allowed to vary. In kinematic 
terms, atoms and bonds are analogous to links and joints of a serial mechanism in which each pair of 
consecutive axes intersects at a link. This implies that the link offsets are zero for all six links which 
allows us to reduce the 6-dimensional problem to a system of 3 polynomials in 3 unknowns. The product 
of all link transformation matrices is the identity matrix, since the end-effector is at the same position 
and orientation as the base link. 

We adopt an approach proposed by D. Parsons [23]. Notation is defined in figure 1. Backbone 
atoms are regarded as points p\, . . . ,p$ G R 3 ; the unknown dihedrals are the angles u)\, . . . ,u>§ about 
axes (p§,pi) and (pi^i,pi) for i = 2, . . . , 6. For readers familiar with the kinematics terminology, the 
Denavit-Hartenberg parameters are 

di = 180° - <pi, di = Li, en = 0,9i = Wj. 

Each of triangles T\ = l\(p\,p2,p§), T2 = A(p2,P3,P4,) and T3 = A{pA,p^>,p%) is fixed for constant 
bond lengths L\,...,L§ and bond angles 01,03,05- Then the lengths of {P21P&), {P2->Pi) and {pi,p%) 
are constant, hence base triangle A(p2,P4,pe) is fixed in space, defining the xy-plane of a coordinate 
frame. Let 9\ be the (dihedral) angle between the plane of A(pi,p2,pe) and the xy-plane. Clearly, for 
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Figure 1: The cyclic molecule. 
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any conformation 0\ is well-defined. Similarly we define angles O2 and #3, as shown in figure 1. We call 
them flap (dihedral) angles to distinguish them from the bond dihedrals. 

Conversely, given lengths Lj, angles 0$ for « = 1, . . . , 6 and flap angles 0j for i = 1, . . . , 3 the coor- 
dinates of all points pi are uniquely determined and hence the bond dihedral angles and the associated 
conformation are all well-defined. We have therefore reduced the problem to computing the three flap 
angles B% which satisfy the constraints on bond angles <j>2, <f>&, 4>q- 

Hence we obtain polynomial system 

a 11 + «12 cos 62 + CC13 cos 6*3 + a 14 cos 62 cos 6*3 + CK15 sin 6*2 sin 6*3 = 0, 
«2i + ^22 cos #3 + «23 cos 9\ + CK24 cos #3 cos 6\ + (X25 sin #3 sin 9\ = 0, 
«3i + «32 cos 6*i + 033 cos 6*2 + 034 cos 6*i cos 6*2 + 035 sin ^i sin ^2 = 0, (12) 

cos 2 6>i + sin 2 9 1 -l = 0, 
cos 2 2 + sin 2 2 -l = 0, 
cos 2 3 + sin 2 3 -l = 0, 

where the are input coefficients. 

This system has Bezout bound of 64 and mixed volume 16; the mixed volume is the exact number of 
complex roots generically as we shall prove below by demonstrating an instance with 16 real roots. Notice 
that 16 is also the exact number of solutions, generically, to the general inverse kinematics problem with 
6 rotational joints (6R). 

For our resultant solver we prefer an equivalent formulation with a smaller number of polynomials, 
obtained by applying the standard transformation to half-angles that gives rational equations in the new 
unknowns tf. 

Q. j2 ^f- 

U = tan ^ : cos B { = — — | , sin 0j = — ^ , i = 1, 2, 3. 

This transformation captures automatically the last three equations in (12). By multiplying both sides of 
the i-th equation by (1 + tj)(l + i|), where k) is a permutation in 5(1, 2, 3), the polynomial system 
becomes 

/i = /?n + /?i2tl + /3i 3 i§ + /3i4*li| + /3i 5 «2i 3 = 

h = 021 + #22*3 + #23*1 + foA^l + 025*3*1 = (13) 
h = 031 + 032*? + 033*2 + 034*1*2 + 035*1*2 = 

where 0jj are input coefficients. The new system has again Bezout bound of 64 and mixed volume 16. 
7.2 Applying the Resultant Solver 

The first instance is a synthetic example for which we fix one feasible conformation with all flap angles 
equal to 90°. All polynomials are multiplied by 8 in order for the coefficients to be all integers, then 0^ 
is the (i,j)-th entry of matrix 

" -9 -1 -1 3 8 " 
-9 -1-13 8 
-9 -1-13 8 

The symmetry of the problem is bound to produce root coordinates of high multiplicity, so we decide to 
follow the first approach to solving the system (sect. 3) and add polynomial 

fo = u + 31*i - 41t 2 + 6U3 
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with randomly selected coefficients. In this system, the 3-fold mixed volumes are 12, 12, 12, 16 hence 
the sparse resultant has total degree 52 and degree 16 in Jq. The resultant matrix is regular and has 
dimension 86, with 30 rows corresponding to / . This is the offline phase; the rest corresponds to the 
online execution of the solver. 

The entire 56 x 56 upper left submatrix is decomposed and is relatively well-conditioned. In the 
30 x 30 matrix polynomial, the leading matrix coefficient is singular within machine precision; two random 
transformations are used but fail to improve significantly the conditioning of the matrix. Therefore the 
generalized eigenproblem routine is called on the 30 x 30 pencil and produces 12 complex solutions, 3 
infinite real solutions and 15 finite real roots. The absolute value of the four polynomials on the candidate 
values lies in [0.6 • 10 _9 ,0.3 • 10~ 3 ] for values that approximate true solutions and in [7.0,3.0 • 10 20 ] for 
spurious answers. Our program computes the true roots to at least 5 digits as seen by comparing with 
the exact solutions computed by Maple V using Grobner bases over the rationals. The latter are 

±(1,1,1), ±(5, -1,-1), ±(-1,5,-1), ±(-1,-1, 5). 

The average CPU time of the online phase on the Sun Sparc 20 of table 2 is 0.4 seconds. 

Usually noise enters in the process that produces the coefficients; this example models this phe- 
nomenon. We consider the cyclohexane molecule which has 6 carbon atoms at equal distances and equal 
bond angles. Starting with the pure cyclohexane, we randomly perturb them by about 10% to obtain 
as the entries of matrix 

-310 959 774 1313 1389 
-365 755 917 1269 1451 . 
-413 837 838 1352 1655 

We used the second approach to define an overconstrained system, namely by hiding variable £3 in the 
coefficient field (sect. 4). The resultant matrix has dimension 16 and is quadratic in £3, whereas the 2-fold 
mixed volumes are all 4 and the sparse resultant has degree 4 + 4 + 4 = 12. 

The monic quadratic polynomial reduces to a 32 x 32 companion matrix on which the standard 
eigendecomposition is applied. After rejecting false candidates, the recovered roots cause the maximum 
absolute value of the input polynomials to be 10~ 5 . We check the computed solutions against those 
obtained by a Grobner bases computation over the integers and observe that each contains at least 8 
correct digits. The total CPU time on a Sun Sparc 20 is 0.2 seconds on average for the online phase. 

Lastly we report on an instance where the input parameters are sufficiently generic to produce 16 real 
roots. The coefficients are given by matrix 

" -13 -1 -1 -1 24 " 
-13 -1 -1 -1 24 . 
-13 -1 -1 -1 24 

We hide £3 and arrive at a resultant matrix of dimension 16, whereas the sparse resultant has degree 12. 
The monic polynomial and the companion matrix are of dimension 32. There are 16 real roots. Four 
of them correspond to eigenvalues of unit geometric multiplicity, while the rest form four groups, each 
corresponding to a triple eigenvalue. For the latter the eigenvectors give us no valid information, so we 
recover the values of t\,t2 by looking at the other solutions and by relying on symmetry arguments. The 
computed roots are correct to at least 7 decimal digits. The average CPU time of the online part is 0.2 
seconds on a Sun Sparc 20. 
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8 Conclusion 



We have examined several computational aspects of sparse elimination theory and, in particular, the 
use of sparse resultant matrices for reducing root-finding to an eigenproblem. A general solver has been 
implemented based on this approach and has been applied successfully to fundamental problems in vision, 
robot kinematics and structural biology. These problems are of moderate size and exhibit sparse structure 
as modeled by the Newton polytopes and the mixed volume. The efficiency and accuracy of our solver 
imply that sparse elimination may be the method of choice for such systems. 

Automating the different ways to deal with numerically unstable inputs will improve the implemen- 
tation. For instance, clustering neighboring eigenvalues and computing the error on the average value 
significantly improves accuracy. A question of practical as well as theoretical interest is to handle the 
case of repeated roots efficiently. 
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